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We present a new algorithm for computing the Lyapunov exponents spectrum based on a matrix 
differential equation. The approach belongs to the so called continuous type, where the rate of 
£SJ ' expansion of perturbations is obtained for all times, and the exponents are reached as the limit at 

infinity. It does not involve exponentially divergent quantities so there is no need of rescaling or 
realigning of the solution. We show the algorithm's advantages and drawbacks using mainly the 
example of a particle moving between two contracting walls. 
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I. INTRODUCTION 



Lyapunov characteristic exponents (LCE) measure the rate of exponential divergence between neighbouring trajec- 
tories in the phase space. The standard method of calculation of LCE for dynamical systems is based on the variational 
equations of the system. However, solving these equations is very difficult or impossible so the determination of LCE 
also needs to be carried out numerically rather than analytically. 

The most popular methods which are used as an effective numerical tool to calculate the Lyapunov spectrum for 
smooth systems relies on periodic Gram-Schmidt orthonormalisation of Lyapunov vectors (solutions of the variational 
equation) to avoid misalignment of all the vectors along the direction of maximal expansion ([l|, Q). 

In some approaches, usually involving a new differential equation instead of the variational one, the procedure of 
re-orthonormalisation is not used [3[ - these are usually called the continuous methods. They are usually found to 
(*C) , be slower than the standard ones due to the underlying equation being more complex than the variational one. A 
comparison of various methods with and without orthogonalisation can be found in 0, [1| and a recent general review 

The main goal of this paper is to present a new algorithm for obtaining the LCE spectrum without the rescaling 
and realigning. This application is a consequence of the equation satisfied by the Lyapunov matrix or operator (see 
below) which was discovered in one of the authors' PhD thesis Q. The particular numerical technique introduced 
here is the first attempt and is open to further development so it still bears the disadvantages of the usual continuous 
• <-h , methods. However, in our opinion, the main advantages of the approach lie in its founding equations and are as 
follows: 

& '■ 

• The whole description of the LCE is embedded in differential geometry from the very beginning, so that and it 
is straightforward to assign any metric to the phase space including one with non-trivial curvature. 

• As the rate of growth is described by an operator (endomorphism) on the tangents space, and the equation it 
satisfies is readily expressed with the absolute derivative, the approach is explicitly covariant. (The exponents 
are obviously invariants then, although their transformation properties still seem to be a live issue, see e.g. [|[.) 

• There is no need for rescaling and realignment, as the main matrix is at most linear with time, and encodes the 
full spectrum of LCE. 

• Since we make no assumptions on the eigenvalues, there are no problems with the degenerate case encountered 
in some other methods. 
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• We rely on a single coordinate-free matrix equation, which reduces the method's overall complexity. 

• The fundamental equation is not an approximation but rather the differential equation satisfied by the so-called 
time-dependent Lyapunov exponents. This opens potential way to analytic studies of the exponents. 

It should be noted that the last points imply a hidden cost (in the current implementation) of diagonalising instead 
of reorthonormalising, due to the complex matrix functions involved. Fortunately, this procedure needs to be carried 
out on symmetric matrices for which it is stable. 

The natural domain of application of this method might be the General Relativity and dynamical systems of 
cosmological origin - already formulated in differential geometric language [914TT1]. Of course this still requires the 
resolution of the question of the time parameter, and natural metric in the whole phase space (not just the configuration 
space which corresponds to the physical space-time). Regardless of that choice, however, the fundamental equation 
of our method will remain the same - whether one chooses to consider the proper interval as the time parameter, 
or find some external time for an eight dimensional phase space associated to the four dimensional space-time. This 
stems from the fact that our approach works on any manifold. 

Here, we wish to focus on the numerical aspect of the method, providing the rough first estimates of its effectiveness. 
This is a natural question, after the theoretical motivation for a given method has been established, namely how well 
it performs numerically. There are obviously many ways of translating the method into code, and we hope for future 
improvements, nevertheless, the presented implementation can be considered a complete, ready-to-use tool. In the 
next sections we review the derivation of the main equation and then proceed to the simple mechanical examples for 
testing and results. 

II. THE BASE SYSTEM OF EQUATIONS 

For a given system of n ordinary differential equations 

x = f{x,t), x£R n , (I) 
the variational equations along a particular solution tp(t) are defined as 

(2) 

Ud - x= V (t) 

and the largest Lyapunov exponent can be defined as 



ox 



»_ . lim MS, (3) 

t— »oo t 



for almost all initial conditions of Z. From now on wc take the norm to be 



\\Z\\:=VZTZ, (4) 

where Z is treated as a column vector, and T denotes transposition. That is to say the metric in the tangent space 
is Euclidean, as is usually assumed for a given physical systems. This needs not be the case, and a fully covariant 
derivation of the main equation can be found in Q. 

The above definitions are intuitively based on the fact that for a constant A, the solution of ([2]) is of the form 

Z = exp(At)Z(0), (5) 

and A max is the greatest real part of the eigenvalues of A. In the simplest case of a symmetric A, the largest exponent 
is exactly the largest eigenvalue. To extend this to the whole spectrum, we note that any solution of ([2} is given in 
terms of the fundamental matrix F so that 

Z = FZ(0), F = AF, F(0) = 1. (6) 

Then (if the limit exists), the exponents are 

/ \n(F T F) \ 

{Ai,...,A„} =spec I lim — . (7) 

\t->oo Zt I 
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Since F T F is a symmetric matrix with non-negative eigenvalues, the logarithm is well defined. The additional factor 
of 2 in the denominator results from the square root in the definition of the norm above. 

As we expect F to diverge exponentially, there is no point in integrating the variational equation in itself, but 
rather to look at the logarithm. To this end we introduce the two matrices M and L: 

M := FF T = F(F T F)F~ 1 =: exp(2L). (8) 

Clearly M has the same eigenvalues as F T F to which it is connected by a similarity transformation, and the eigenvalues 
of L behave as tX for large times 

{Ai, . . . , A„} = spec ( lim ^ ) . (9) 

V t^oo t J 

That is why we call L the Lyapunov matrix. 

To derive the differential equation satisfied by L we start with the derivative of M and use the property of the 
matrix (operator) exponential 

l 

M = AFF T + FF T A T = e 2L / e- 2aL 2Le 2aL da 

(10) 



o 

i 



- 2L Ae 2L + A T = 2 j e- 2a ^Lda 



where we have introduced a concise notation for the the adjoint of L acting on any matrix X as 

[L(X) = [L,X], (11) 

and used its property 

e~ L Xe L = e~ [L X. (12) 
Next, the integral is evaluated taking the integrand as a formal power series in [L 

e-^A + A T = 2^^L, (13) 

where the fraction is understood as a power series also, so that there are in fact no negative powers of [L. Alternatively 
one could justify the above by stating that the function 

1 - e~ 2u 

Mu) = (14) 
u 

is well behaved on the spectrum of [L which is contained in R. As ipi is never zero for a real argument, we can invert 
the operator on the left-hand side of flT3"|) to get 



[L 



e-^ L (9 + uj) + 9-Lj) =L, (15) 



1 - e- 2 ^ 

where the symmetric and antisymmetric parts of A are 

6 = \{A + A T ), u = \{A T -A T ). (16) 

This allows for the final simplification to 

L = [L coth ([£) 9 — [L, uj] (17) 

The function ip%{u) = ucoth(u) should be understood as the appropriate limit at u = 0, so that it is well behaved for 
all real arguments. As was proven in [7J, the above equation is essentially the same in general coordinates: 

V v L = MIL)0-[L,lu], (18) 
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where v is the vector field associated with ([T]) , and V is the covariant derivative. 

Note that in this form it is especially easy to obtain the known result for the sum of the exponents. Since trace of 
any commutator is zero, the only term left is the "constant" term of ip2 which is 1 (or rather the identity operator) 
so that 



where V is the volume of the parallelopiped formed by n independent variation vectors. 

Another simple consequence occurs when the 9 matrix is zero, the whole equation becomes a Lax equation 



which preserves the spectrum over time, so that L/t tends to zero at infinity. Another way of looking at it is that it 
is a linear equation in L and the matrix of coefficients [lo is antisymmetric in the adjoint representation, so that the 
evolution is orthogonal and the matrix norm (Frobenius norm to be exact) of L is constant which means L/t tends 
to the zero matrix. The simplest example of this is the harmonic oscillator or any critical point of the centre type. 
The variations are then vectors of constant length and the evolution becomes a pure rotation. The authors are not 
aware of any complex or non-linear system that would exhibit such simple behaviour. Already for the mathematical 
pendulum such picture is achieved only asymptotically for solutions around its stable critical point. One could expect 
that a system with identically zero exponents might not be "interesting" enough to incur this kind of research. 

We have thus arrived at a dynamical system determining L, with right-hand side being given as operations of the 
adjoint of L on time dependent (through the particular solution) matrices 8 and uj. The next section deals with the 
practical application of the above equation. 



The main difficulty in using (|17p is the evaluation of the function of the adjoint operator. Since we will be integrating 
the equation to obtain the elements of the matrix L, it would be best to have the right-hand side as an explicit 
expression in those elements. This can be done for the n = 2 case, but already for n = 3 one has dozens of terms on 
the right, and for higher dimensions the number of terms is simply too large for such an approach to be of practical 
value. An alternative (although equivalent) dynamical system formulation for the mentioned low dimensionalities 
have been studied in [3j, but again the complexity of the equations increases so fast with the dimension that the 
practical value is questionable. Our method, on the other hand, can be made to rely on the same equation for all 
dimensions, and the only complexity encountered will be the diagonalisation of a symmetric matrix of increasing size. 

Another problem lies in the properties of the ip2 (u) function which, although finite for real arguments, has poles at 
u = ±iir. This means that a series approximation is useless, as it would converge only for eigenvalues smaller than n 
in absolute value, whereas we expect them to grow linearly with time and need the results for t — > oo. On the other 
hand, ip2{u) — > \u\ for large u but, unfortunately, the adjoint operator always has n eigenvalues equal to zero, and for 
Hamiltonian systems it is also expected that two eigenvalues tend to zero. 

Thus, as we require the knowledge of tp2 for virtually any symmetric matrix L, and we are going to integrate the 
equation numerically anyway, we will resort to numerical method for this problem. Because the matrix L is symmetric 
(Hcrmitian in an appropriate setting) so is its adjoint [L, and the best numerical procedure to evaluate its functions 
is by direct diagonalisation 12]. Obviously this is the main disadvantage of the implementation method as even for 
symmetric matrices, finding the eigenvalues and all the eigenvectors is time-consuming. So far the authors have only 
been able to find one alternative routine which is to numerically integrate not L itself but rather the diagonal matrix 
of its eigenvalues and the accompanying transformation matrix of eigenvectors. However, due to the increased number 
of matrix multiplications the latter method does not seem any faster than the former. 

With this in mind, let us see how the diagonal form of L simplifies the equation. First, we need to regard [L as an 
operator, and since it is acting on matrices we will adopt a representation where any n x n matrix becomes a n 2 x 1 
matrix, i.e. a n 2 element vector constructed by writing all the elements of successive rows as one column. [L is then 
a n 2 x n 2 matrix. Fortunately, one does not need to diagonalise [L but only L itself. As can be found by direct 
calculation, the eigenvalues of the adjoin are all the differences of the eigenvalues of L. For example 




(19) 



L = [w,L], 



(20) 



III. EXEMPLARY IMPLEMENTATION 




(21) 
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where the subscript D denotes "diagonal" . 

Now let us assume we also have the transformation matrix such that 

L = QL D Q T . 



(22) 



Then, instead of bringing the whole equation to the eigenbasis of L, one can only deal with the 9 matrix in the 
following way 



M[L)0 = QM[Ld)Q T = QM[Ld)0Q T , := Q T 9Q. 



(23) 



Of course, the other term of equation ([T71) can be evaluated as the standard commutator. For the above example the 
ip2 part would be 



M[Ld) = 
and converting 9 to a vector we get 



a o o 

ML1-L2) 






1 

M^2- 



(24) 



Li) 



hi 622 



= M[Li 



(9n\ 




( - \ 


S\2 




^2(Ai 2 )6»i2 


921 




^ 2 (A 21 )0 21 


\9~22) 




V O22 J 



(25) 



which corresponds to the usual 2 by 2 matrix of 



9ll _ M^J2)9l2 
■0 2 (A 2 l)021 9~22 



(26) 



where Ay := Li — Lj are the differences of the eigenvalues of L. In general the appropriate (n x n) matrix elements 
read 



(27) 



and this matrix needs to be transformed back to the original basis according to ()23[) before being used in the main 
equation. 

The matrix elements of L will, in general, grow linearly with time. This is of course a huge reduction when compared 
with the exponential growth of the perturbations, but one might want to make them behave even better by taking 
time into account with 



L = (t + l)A. 



(28) 



The reason for the additional 1 is that we will specify the initial conditions at t = and we want to avoid dividing 
by zero in the numeric procedure and the limit at infinity (if it exists) is not affected by this change. Of course, in 
the case of the autonomous systems any value of t can be chosen as initial (at the level of the particular solution) but 
the non-autonomous case might require a particular value, which can be dealt with in a similar manner. 
We now have 



A = 



f 



t + 1 



(M(t + l)[A)0-A)-[A,w] 



(29) 



As for the initial conditions, the fundamental matrix F is equal to the identity matrix at t = 0, so that L(0) = 
which in turn implies A(0) = 0. 

We are now ready to state the general steps of the proposed implementation. Choosing a specific numerical routine 
to obtain the solution for a time step h at each t these are 

1. Obtain the particular solution ip(t), calculate the Jacobian matrix A at (p and, from it, the two matrices 9(t) 
and oj(t). Start with A(0) = 



2. Find the eigenvalues A, and eigenvectors Q of A 
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3. Transform 6(t) to 6 = Q T 6Q. 

4. Compute the auxiliary matrix /Jy = [ip((t + l)[Aa)0]ij = tp2((t + l)Aij)9ij 

5. Take the derivative to be A(t) — (Q f3Q T — A)/(t + 1) — [A,uj], and use it to integrate the solution to the next 
step A(t + h). 

6. Repeat steps 2-5 until large enough time is reached. 

7. The "time-dependent" Lyapunov exponents at time t are simply (1 + j)Aj. 

The relation in the last point stems from the rescaling (1281) and, as can be seen, is only important for small values of 
t. 



IV. EXAMPLES AND COMPARISON 



As suggested in the preceding section, one expects this method of obtaining the Lyapunov spectrum to be relatively 
slow. In order to sec that, we decided to compare it with the standard algorithm based on direct integration of 
the variational equation and Gram- Schmidt rescaling at each step [![. Although usually the rescaling is required 
after times of order 1, we particularly wish to study a system with increasing speed of oscillations and frequent 
renormalisation will become a necessity. In other words, we want to give the standard method a "head start" when 
it comes to precision. 

For numerical integration we chose the modified midpoint method (with 4 divisions of the whole timestep H = 4/i) , 
which has the advantage of evaluating the derivative fewer times than the standard Runge-Kutta routine with the 
same accuracy. Since we intended a simple comparison on equal footing, we did not try to optimise either of the 
algorithms and wrote the whole code in Wolfram's Mathematica due to the ease of manipulation of the involved 
quantities and operations (e.g. matrices and outer products). 

The most straightforward comparison is for the simplest, i.e. linear dynamical system, for which the main equation 
([T]) is the same as the variational one @ , with constant matrix A. The exponents are then known to be the real parts 
of the eigenvalues of A. To include all kinds of behaviour, we took a matrix which has a block form 



A = 



(l -8 0\ 

12 1 

| 2 

V° 1 y 



(30) 



whose eigenvalues are {1 + 4i\/6, 1 — 4i-\/6, \ + V2, \ — V2}, so that the LCE are approximately equal to 
{1,1,1.91421,-0.914214}. 

The basic timestep was taken to be h = 0.01 and the time to run from to 1000. The results are depicted in figures [1] 
and[2]with the horizontal axis representing the inverse time 1 jt so that the sought for limit at infinity becomes the value 
at zero which is often clearly seen from the trend of the curves. We note that our method required 59 seconds, whereas 
the standard one only took 17 seconds. The final values of LCE (t = 1000) were {1.00010, 0.99990, 1.91427, -0.91427} 
and {0.99933,0.99970, 1.91372, —0.91372}, respectively. The shape of the curves is different, which is to be expected 
because the matrix A measures the true growth of the variation vectors at each point of time, and the other method 
provides more and more accurate approximations to the limit values of the spectrum. 

Before we go on to the central example, which explores a system with accelerating oscillations, we will present 
numerical results for the system which is synonymous with chaos, namely the Lorenz system (l3| . 

The equations read 

x = <j(y — x) 

y = x(p-z)- y (31) 
z — xy — j3z, 

where we took er = 10, j3 = 8/3 and p = 28, and integrated the equation for the initial conditions of x — y = 
z = 10 from t = to 1000. Next we integrated the respective methods for the exponents with the timestep of 
0.001. Our method took about 591 seconds and the result is shown in figure [3] with the final value of the spectrum 
{-14.5676, -0.0010078, 0.901931}. Note that we have shifted the lowest exponent by +15, so that all three could be 
presented on the same plot with enough detail. The standard method took about 152 seconds and its outcome is 
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Figure 1: The Lyapunov exponents for the linear system obtained with the presented method. 
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Figure 2: The Lyapunov exponents for the linear system obtained with the standard method. 



shown in figure |4] with the final values of {—14.5664,0.00173051,0.897989} (we have shifted the graph in the same 
manner as before). 

One could note that there is less overall variation of the time-dependent exponents in our method similarly to the 
previous example. A good estimate of precision would be to calculate the sum of the exponents, which, in this system, 
should be exactly equal to -41/3. The difference between that and the numerical estimates were: for the standard 
method (at t = 1000) -9.04 x 10" 7 , and for our method —9.67 x 10 11 - a much better result. This seems to be the 
usual picture for the continuous methods which trade computing time for precision. 

A presentation of some more complicated, including both integrable and chaotic, examples can be found in Q, and 
for such systems also, we observe the concordance of final results and the speed discrepancy. In order to see how one 
could benefit from the new method we have to turn to another class of models, ones for which the exponents present 
oscillatory behaviour. We found that for artificial systems with accelerated oscillations our method performs better 
and present here a simple physical model which exhibits such property. 

Consider a ball moving between two walls which are moving towards each other and assume that the ball bounces 
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Figure 3: The Lyapunov exponents for the Lorenz system obtained with the new method. The middle exponent has been 
shifted by +15 to show the details with a better scale. 
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Figure 4: The Lyapunov exponents for the Lorenz system obtained with the standard method. The middle exponent has been 
shifted by +15 to show the details with a better scale. 



off with perfect elasticity. As the distance between the walls decreases and the speed of the ball increases it takes less 
and less time for each bounce cycle. In order to model this analytically, without resorting to infinite square potential 
well we will take the following Hamiltonian system 

H =\p 2 + U {j(r))-> U(x) = ^Tt & nh{x)\ (32) 

which depends on time explicitly via the function / whose meaning is seen as follows: the area hyperbolic tangent is 
infinite for x = 1, so that the potential becomes infinite for the position variable q = /(£), so that / is simply half the 
distance between the walls. In particular we take it to be 

no - m 
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so that it decreases from 1 to \ slowing down but never stopping. The reason for this is that we want the system not 
to end in a finite time, and also that the worse behaved / is the faster the numerical integration of the main system 
itself will fail. The initial shape of the potential and the systems setup is depicted in figure [5] 
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Figure 5: The potential modelling two contracting walls. The vertical lines correspond to the impenetrable barriers of "infinite" 
U. (The vertical position of the ball has no meaning.) 

The slowing down of the walls, and the particular shape of U allows us to find rigorous bounds on the Lyapunov 
exponents. First we note, that the vector tangent to a trajectory in the phase space, i.e. a vector whose components 
are simply the components of F from |T]), is always (for any dynamical system) a solution of the variational equation. 
That means that just by measuring its length we can estimate the largest exponent, since for almost all initial 
conditions the resulting evolution is dominated by the largest exponent. Second, the system is Hamiltonian and it 
must have two exponents of the same magnitude but different signs, so analysing this particular vector will give us 
all the information regardless of the initial condition. Thus we have to find the following quantity 



A( t ) = -ln, 



P 



(34) 



as a function of time, which boils down to finding bounds on the velocity and acceleration of the bouncing ball. 
As the Hamiltonian depends on time explicitly, the energy is not conserved, but instead we have 



(35) 



The velocity has its local maxima at q = when all the energy is in the kinetic term, and we are lead to define a 
virtual maximal velocity by equating the energy at any given time to a kinetic term 



f 



-v M := E, 



(36) 



we will take the positive sign of vm, and assume it is non-decreasing as the physical setup suggest. Differentiating 
the above we get 



vmvm = -U' 



q- 



f 



J J V 2 ' 

Let us go back to the equation of motion for the momentum variable which reads 

' — G)7- 



(37) 



(38) 
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and substitute that into the previous equation to get 

/ 

VMVM=I>qj- (39) 

As mentioned above / decreases very slowly at late times, which is when we estimate the exponents anyway. We thus 
assume, that / is small enough for the fraction on the right-hand side to be considered constant over one cycle - that 
is over the time T in which the ball moves from the centre up the potential wall and back to the centre. This time 
will get shorter and shorter, but also / will get closer and closer to zero. The standard problem of the elastic ball 
and infinitely hard walls shows that the speed transfer at each bounce is of the order of / so the (virtual) maximal 
velocity will change as slowly as / and we are entitled to average the equations over one cycle: 

. t+T . / t+T \ 

vmvm = \- f J pq dt = [qp]l +T - J P 2 dt < - f -v 2 Ml (40) 

where we integrated by parts and used the fact that at the beginning and end of the cycle q = 0, and also that the 
momentum is never greater than the maximal velocity. 

We are thus left with a bound in the form of a differential equation, and since all the quantities involved are positive 
and non-decreasing (— / = |/|), its solution will be the bound for the velocities 

q < v M < -j-, (41) 

with the "initial" value of vm taken at sufficiently large t so that / is small. 

Similar considerations can be carried out for the acceleration p, only this time we have to introduce a virtual points 
of return g^j, that is the point at which all the energy is in the potential term 

U (^j := E, (42) 

since at the real turning points the acceleration reaches its local maxima. Note that this is not the same as vm which 
describes the slow growth of the consecutive maxima of q and not the maxima of its slope. 
This definition allows us to express qn as a function of vm (via energy) 



q r = ftanh(v M ), (43) 



and the acceleration is 



P < O-M 



-U' 



/artanh 

P~q 2 R ~ J(l-tanh 2 (w A f)) (44) 



VM 



VMO ,2 ( V M0 

dM = -p- cosh \ 

Bringing the two results together we see that 

q 2 +P 2 <v 2 M + a\ { < oo, (45) 



because the function / does not tend to zero, and by (|34[) both Lyapunov exponents must in turn be zero themselves. 
We also recognise that the hyperbolic cosine factor in ojf could produce a nonzero exponents if / were to decrease to 
zero as 1/t. 

Let us now turn to the numerical results of both methods for this system. As initial conditions we take q(0) = 
and p(0) = 1. For the same time step h — 0.01 we see in figure [6] that the new method predicts the values correctly, 
integrating for 59 seconds from t = to t = 1000. However for the standard one, as shown in figure [71 we see the 
exponents diverging from zero, for the same time limits, and integrating time of 17 seconds. It turns out h = 0.005 
also gives divergent results and the correct behaviour is recovered for h = 0.001, shown in figure for which the 
routine takes 171 seconds. 
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Figure 6: The Lyapunov spectrum obtained with the new method for two contracting walls (timestep h — 0.01). 
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Figure 7: The Lyapunov spectrum obtained with the standard method for two contracting walls (time step h — 0.01). 

V. CONCLUSION 

We have presented a new algorithm for evaluation of the Lyapunov spectrum, emerging in the context of differential 
geometric description of complex dynamical systems. This description seems especially suitable for systems found in 
General Relativity like, e.g., chaotic geodesic motion [141 ], The main advantage of the base method is its covariant 
nature and concise, albeit explicit, matrix equations that promise more analytic results in the future. Also, this allows 
for study of curved phase spaces and general dynamical systems - not only autonomous or Hamiltonian ones. 

The main differential equation, can be numerically integrated, giving a simple immediate algorithm for the com- 
putation of the Lyapunov characteristic exponents. It is in general slower than the standard algorithm (based on 
Gram-Schmidt orthogonalisation), but the first numerical test suggest it works betters in systems with increasing 
frequency of (pseudo-)oscillations. We show this on the example of a simple mechanical system - a ball bouncing 
between two contracting walls. 

Although in low dimensions the main equation can be cast into an explicit form (with respect to the unknown 
variables), in general the numerical integration requires diagonalisation at each step, which is the main disadvantage 
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Figure 8: The Lyapunov spectrum obtained with the standard method for two contracting walls (time step h — 0.001). 

of the method and the reason of its low speed. We hope to present a more developed algorithm without this problem 
in the future. 
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